classdef Gauss < FilterFather
    properties
    end
    methods
        
        function this = Gauss(guiData)
            this.calculateFilter(guiData);
%             minEdge = min([guiData.Fp 0.1])*2*pi;
%             maxEdge = max([guiData.Fp guiData.Fa])*2*pi;
%             this.calcData([minEdge maxEdge]);
        end

        function calculateFilter(obj,guiData)
%             Aa = guiData.Aa;
%             Ap = guiData.Ap;
            obj.FilterType = guiData.filterType;
            guiData.Fa = 10*guiData.Fp;
            fixOrder = guiData.order;
            minEdge = min([guiData.Fp guiData.Fa])*2*pi;
            maxEdge = max([guiData.Fp guiData.Fa])*2*pi;
%             Wgd = guiData.Fp*2*pi;
%             Wa = guiData.Fa*2*pi;
%             DesnormValue = guiData.DenormVal;
%             DesnormType = guiData.DenormType;
%             fixOrder = guiData.fixOrder;
%             MaxOrder = guiData.order;
%             MaxQ = guiData.maxQ;

            tau =  guiData.Delay;
            tol =  guiData.Tolerance/100;
           
            
           % plantilla de default
%             Ap=3;
%             Aa=25;
%             Wa=1.5;
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
%                 N = obj.GaussOrder(Wgd,Ap,Aa,tau,tol);
%                 
%                 if MaxOrder == 0
%                     obj.Order = N;
%                 else 
%                     if(N > MaxOrder || fixOrder == true)
%                         N = MaxOrder;
%                     end
%                     obj.Order = N;
%                 end
%                 
                %Gauss normalizado
%                 [z p k] = ZPKadjust(z,p,0);
    
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            start = max([1 fixOrder]);
            works = 0;
            for N = start: 15
                [z p k] = obj.GaussZPK(N);
                obj.K = abs(real(prod(obj.Poles)));
                num = Poly(z);
                den = Poly(p);
                obj.TransferFunction = tf(num,den);
                obj.calcData([1 2]);
                Wgd = obj.GroupDelay(1)/tau;
                [num,den] = lp2lp(num,den,Wgd);
                obj.Zeros = roots(num);
                obj.Poles = roots(den);
                w = 2*pi*linspace(0.1,1,50000);
                H = tf(num,den);
                [Mag] = bode(H,w);
                H = H/(max(Mag));
                obj.TransferFunction = H;
                obj.TransferFunction = obj.TransferFunction*10^(guiData.gain/20);
                obj.calcData([minEdge maxEdge]);
%                 obj.calcData ([1 2]);
                k = find(obj.Frec > guiData.Fp, 1, 'first');
                if (obj.GroupDelay(k-1) > tau - tau*tol) || fixOrder
                    works = 1;
                    break
                end
            end
            
            if ~works
                obj.Order = 40;
            else
                obj.Order = N;
            end

            

            
        end
        
        function N = GaussOrder( obj, Wgd, Ap, Aa, tau, tol)
        %Wn: frecuencia de corte en el retardo de grupo 
        %tau: retardo de grupo
        %tol: tolerancia para el retardo de grupo
            MaxOrden = 15;    % numero masonico, edad en la que murio cristo, etc..
            N = MaxOrden;
            
                for N = 1:1;
                
                    [z p k] = obj.GaussZPK(N, Ap);
                    [z p k] = ZPKadjust(z, p, 0);

                    H = tf(poly(z),poly(p));
                    W = logspace(log10(Wgd/1000), log10(Wgd*100), 8000);
                    [Mag,Phase] = bode(H,W);
                    Phase = deg2rad(squeeze(Phase));
                    dW = diff(W);
                    dph = diff(Phase)';
                    GD = -dph./dW;
                    W = W(2:end);
    
                    tau0 = GD(1);   %GD inicial, dividiendo por el mismi, normalizo en 1
                    t = find(W > Wgd*tau, 1, 'first');
                    if ((GD(t)/tau0) >= tol)
                       % figure(2)       %DEBUGGING
                       % semilogx(W/tau,tau*GD/tau0)
                        break;
                    end
                
                end
                if(N==MaxOrden)
                    errordlg('Maximum order for Gauss approximation reached. (n = 15)'); 
                end
            
        end
        
        
        function [z p k] = GaussZPK(obj,n)

        %|H(jw)|^2 = 1/(1 + sum(k=1:n){ (log(ro)^k / k!) * (w/wp)^2k})
        %ro = 10^(-Ap/10);
        
%             ro = 10^(-Ap/10);
            gamma = 0.995;
            denominador = zeros(1,2*n+1);
%             denominador(1) = 1;
     
            for k = 1:n
                denominador(end - 2*k) = (-1)^k*(gamma^k) / prod(1:k);
            end
%             denominador = fliplr(denominador);
            denominador(end) = denominador(end) + 1;
            p = roots(denominador);
            p = p(real(p)<0);
            z = [];
            k = abs(prod(p));
        end
        
        
    end
    
end

